{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2ec58e18-ef54-45aa-911b-0ee33c8415e3",
   "metadata": {},
   "source": [
    "Chapter 03\n",
    "\n",
    "# Python自定义函数计算矩阵乘法，外积视角\n",
    "《线性代数》 | 鸢尾花书：数学不难"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "991e32be-cde0-4476-b686-73ce804da7e8",
   "metadata": {},
   "source": [
    "这段代码从**矩阵乘法的“第二视角”——列-行外积展开（Outer Product Expansion）**出发，对标准矩阵乘法进行了实现。在第一视角中，我们强调**每个元素是如何由行与列的点积计算得出**；而这里的“第二视角”强调的是：**整个矩阵乘积是多个秩1矩阵（外积）的线性叠加**。下面从数学角度详细解析，并穿插相应的公式表示。\n",
    "\n",
    "---\n",
    "\n",
    "### **1. 数学背景：外积展开视角的矩阵乘法**\n",
    "\n",
    "给定两个矩阵：\n",
    "\n",
    "- $A \\in \\mathbb{R}^{m \\times p}$\n",
    "- $B \\in \\mathbb{R}^{p \\times n}$\n",
    "\n",
    "设：\n",
    "- $A$ 的列向量表示为 $A = \\left[ \\mathbf{a}^{(1)}, \\mathbf{a}^{(2)}, \\dots, \\mathbf{a}^{(p)} \\right]$, 其中 $\\mathbf{a}^{(k)} \\in \\mathbb{R}^m$\n",
    "- $B$ 的行向量表示为 $B = \\begin{bmatrix} (\\mathbf{b}^{(1)})^T \\\\ (\\mathbf{b}^{(2)})^T \\\\ \\vdots \\\\ (\\mathbf{b}^{(p)})^T \\end{bmatrix}$，其中 $\\mathbf{b}^{(k)} \\in \\mathbb{R}^n$\n",
    "\n",
    "那么，矩阵乘积 $C = AB$ 可以表示为**外积求和的形式**：\n",
    "\n",
    "$$\n",
    "C = AB = \\sum_{k=1}^{p} \\mathbf{a}^{(k)} (\\mathbf{b}^{(k)})^T\n",
    "$$\n",
    "\n",
    "其中每一项 $\\mathbf{a}^{(k)} (\\mathbf{b}^{(k)})^T$ 是一个 $m \\times n$ 的**秩1矩阵**（即行 rank 为1 的矩阵），代表了矩阵 $A$ 第 $k$ 列与矩阵 $B$ 第 $k$ 行之间的**外积**。\n",
    "\n",
    "- $\\mathbf{a}^{(k)}$ 是一个列向量 $\\in \\mathbb{R}^{m \\times 1}$\n",
    "- $(\\mathbf{b}^{(k)})^T$ 是一个行向量 $\\in \\mathbb{R}^{1 \\times n}$\n",
    "- 它们的乘积 $\\mathbf{a}^{(k)} (\\mathbf{b}^{(k)})^T \\in \\mathbb{R}^{m \\times n}$\n",
    "\n",
    "这个视角强调了**矩阵乘积是一系列“结构简单”的秩1矩阵的叠加**，在神经网络、图像重构和张量分解中有广泛的应用。\n",
    "\n",
    "---\n",
    "\n",
    "### **2. 代码实现逻辑**\n",
    "\n",
    "- 程序首先检查维度匹配，即确保 $A$ 的列数等于 $B$ 的行数。\n",
    "- 初始化结果矩阵 $C$ 为零矩阵，维度为 $m \\times n$。\n",
    "- 对于每一个 $k = 0, 1, \\dots, p-1$：\n",
    "  - 从矩阵 $A$ 中提取第 $k$ 列向量 $a_k = A[:, [k]]$，形状为 $(m, 1)$\n",
    "  - 从矩阵 $B$ 中提取第 $k$ 行向量 $b_k = B[[k], :]$，形状为 $(1, n)$\n",
    "  - 计算外积：$a_k \\cdot b_k$（用 `@` 运算符表示矩阵乘法）\n",
    "  - 将该秩1矩阵加到 $C$ 中：$C \\mathrel{+}= a_k \\cdot b_k$\n",
    "\n",
    "最终，$C$ 中累加了所有秩1矩阵之和，得到结果 $AB$。\n",
    "\n",
    "---\n",
    "\n",
    "### **3. 实例计算分析**\n",
    "\n",
    "仍然考虑矩阵：\n",
    "\n",
    "$$\n",
    "A = \\begin{bmatrix}\n",
    "1 & 2 & 3 \\\\\n",
    "4 & 5 & 6\n",
    "\\end{bmatrix}, \\quad B = A^T = \\begin{bmatrix}\n",
    "1 & 4 \\\\\n",
    "2 & 5 \\\\\n",
    "3 & 6\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "则乘积 $C = AB$ 是：\n",
    "\n",
    "$$\n",
    "AB = \\sum_{k=1}^{3} \\mathbf{a}^{(k)} (\\mathbf{b}^{(k)})^T\n",
    "$$\n",
    "\n",
    "每一项如下（逐项展示）：\n",
    "\n",
    "- 第1项（$k=1$）：\n",
    "  $$\n",
    "  \\mathbf{a}^{(1)} = \\begin{bmatrix} 1 \\\\ 4 \\end{bmatrix}, \\quad \\mathbf{b}^{(1)} = \\begin{bmatrix} 1 & 4 \\end{bmatrix} \\Rightarrow\n",
    "  \\mathbf{a}^{(1)} (\\mathbf{b}^{(1)})^T =\n",
    "  \\begin{bmatrix}\n",
    "  1 \\cdot 1 & 1 \\cdot 4 \\\\\n",
    "  4 \\cdot 1 & 4 \\cdot 4\n",
    "  \\end{bmatrix}\n",
    "  =\n",
    "  \\begin{bmatrix}\n",
    "  1 & 4 \\\\\n",
    "  4 & 16\n",
    "  \\end{bmatrix}\n",
    "  $$\n",
    "\n",
    "- 第2项（$k=2$）：\n",
    "  $$\n",
    "  \\mathbf{a}^{(2)} = \\begin{bmatrix} 2 \\\\ 5 \\end{bmatrix}, \\quad \\mathbf{b}^{(2)} = \\begin{bmatrix} 2 & 5 \\end{bmatrix} \\Rightarrow\n",
    "  \\mathbf{a}^{(2)} (\\mathbf{b}^{(2)})^T =\n",
    "  \\begin{bmatrix}\n",
    "  4 & 10 \\\\\n",
    "  10 & 25\n",
    "  \\end{bmatrix}\n",
    "  $$\n",
    "\n",
    "- 第3项（$k=3$）：\n",
    "  $$\n",
    "  \\mathbf{a}^{(3)} = \\begin{bmatrix} 3 \\\\ 6 \\end{bmatrix}, \\quad \\mathbf{b}^{(3)} = \\begin{bmatrix} 3 & 6 \\end{bmatrix} \\Rightarrow\n",
    "  \\mathbf{a}^{(3)} (\\mathbf{b}^{(3)})^T =\n",
    "  \\begin{bmatrix}\n",
    "  9 & 18 \\\\\n",
    "  18 & 36\n",
    "  \\end{bmatrix}\n",
    "  $$\n",
    "\n",
    "求和得到最终结果：\n",
    "\n",
    "$$\n",
    "AB = \\begin{bmatrix}\n",
    "1 + 4 + 9 & 4 + 10 + 18 \\\\\n",
    "4 + 10 + 18 & 16 + 25 + 36\n",
    "\\end{bmatrix}\n",
    "= \\begin{bmatrix}\n",
    "14 & 32 \\\\\n",
    "32 & 77\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "这与标准矩阵乘法一致，但构造方式完全不同，强调结构性。\n",
    "\n",
    "---\n",
    "\n",
    "### **4. 总结**\n",
    "\n",
    "这段代码展示了矩阵乘法的**外积展开法（Outer Product Expansion）**的数学实质：\n",
    "\n",
    "- 如果说第一视角强调**元素的点积计算**（局部视角），\n",
    "- 那么第二视角强调的是**矩阵整体是秩1矩阵的叠加**（全局视角）。\n",
    "\n",
    "这种方法在许多应用中极具启发性，特别是在需要对矩阵进行**低秩近似（low-rank approximation）**、**奇异值分解（SVD）**、**神经网络权重分解**等操作时。它揭示了矩阵乘法中“结构的本质”：一个复杂的矩阵结果可以由多个简单的 rank-1 结构组合而成。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88c630a4-9ef3-412a-bd5d-0bbcd241263e",
   "metadata": {},
   "source": [
    "## 初始化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e2f2e5b3-ed89-42fa-ac0c-ca3d5fce50df",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81b93b87-d4d0-4f25-99f5-10eab114133a",
   "metadata": {},
   "source": [
    "## 自定义函数，矩阵乘法第二视角"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "15069f4e-9172-4627-8167-9889980d2c59",
   "metadata": {},
   "outputs": [],
   "source": [
    "def matrix_multiplication_outer(A, B):\n",
    "    \n",
    "    # 获取矩阵 A 和 B 的形状\n",
    "    m, p_A = A.shape\n",
    "    p_B, n = B.shape\n",
    "\n",
    "    # 检测矩阵形状是否符合矩阵乘法规则\n",
    "    if p_A != p_B:\n",
    "        raise ValueError('Dimensions do not match')\n",
    "\n",
    "    # 初始化结果矩阵 C，形状 (m, n)，初始值设为 0\n",
    "    C = np.zeros((m, n))\n",
    "\n",
    "    for k in range(p_A):\n",
    "        a_k = A[:, [k]]      # A 的第 k 列，二维数组\n",
    "        b_k = B[[k], :]      # B 的第 k 行，二维数组\n",
    "        C += (a_k @ b_k)     # 每次计算外积并加到结果中\n",
    "\n",
    "    return C"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e0b13c6-2d63-4b48-b980-6131310226c8",
   "metadata": {},
   "source": [
    "## 矩阵乘法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3ad79134-302f-4dce-a967-f1877f8360fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2, 3],\n",
       "       [4, 5, 6]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[1, 2, 3],\n",
    "              [4, 5, 6]])\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6151fbf4-afbe-48f1-82f9-e82ff35c2c93",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 4],\n",
       "       [2, 5],\n",
       "       [3, 6]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = A.T\n",
    "B"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61315e03-1f9f-4457-abbd-b827c4beed84",
   "metadata": {},
   "source": [
    "## 矩阵乘法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "40170372-8d49-4525-9b3a-c0c115bef981",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[14., 32.],\n",
       "       [32., 77.]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "matrix_multiplication_outer(A, B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "8adb8924-ee53-4204-af76-e26fbc024feb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[17., 22., 27.],\n",
       "       [22., 29., 36.],\n",
       "       [27., 36., 45.]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "matrix_multiplication_outer(B, A)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25467aae-3221-4bb9-ae3b-d0a86b7f03c2",
   "metadata": {},
   "source": [
    "作者\t**生姜DrGinger**  \n",
    "脚本\t**生姜DrGinger**  \n",
    "视频\t**崔崔CuiCui**  \n",
    "开源资源\t[**GitHub**](https://github.com/Visualize-ML)  \n",
    "平台\t[**油管**](https://www.youtube.com/@DrGinger_Jiang)\t\t\n",
    "\t\t[**iris小课堂**](https://space.bilibili.com/3546865719052873)\t\t\n",
    "\t\t[**生姜DrGinger**](https://space.bilibili.com/513194466)  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:base] *",
   "language": "python",
   "name": "conda-base-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
